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Abstract 

A  mathematical  model  has  been  developed  for  a  lead/acid  cell  with  immobilized  electrolyte,  using  the  transport  equations 
for  concentrated  electrolyte  solutions  and  experimentally  determined  kinetic  equations  for  the  two  porous  electrodes.  The 
model  equations  have  been  solved  by  means  of  an  algorithm,  which  gives,  after  discretization  of  the  differential  equations,  a 
single  tridiagonal  coefficient  matrix  for  the  concentration  profile  over  the  whole  multi-region  system  consisting  of  a  positive 
electrode,  a  separator  and  a  negative  electrode,  respectively.  The  model  can  be  used  for  optimization  of  the  actual  type  of 
lead/acid  cell  for  different  applications.  The  numerically  predicted  results  are  in  fair  agreement  with  experimental  data  for 
acid  and  lead  sulfate  concentration  profiles  in  a  cell,  although  the  experimentally  measured  acid  concentration  in  the  negative 
electrode  after  discharge  is  generally  lower  than  the  calculated  value,  and  the  predicted  discharge  time  is  lower  than  the 
experimental  ones  at  high  discharge  rates. 
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1.  Introduction 

Mathematical  modelling  of  the  lead/acid  battery  is 
a  powerful  tool  in  the  attempt  to  analyse,  predict  and 
optimise  its  performance,  both  in  its  traditional  ap¬ 
plications  as  a  car  battery  for  starting,  lighting  and 
ignition  or  as  a  truck  battery,  as  well  as  in  new  ap¬ 
plications  such  as  uninterrupted  power  supply,  load 
levelling  and  electric  cars.  Mathematical  models  for 
the  lead/acid  cell  have  been  reported  earlier.  A  porous 
electrode  model  for  the  positive  plate  was  extended 
first  to  the  negative  plate  and  then  to  the  complete 
cell  by  Micka  and  Rousar  [1],  They  assumed  that  the 
electrolyte  between  the  two  plates  is  well  mixed  by 
free  convection  and  they  calculated  its  mean  concen¬ 
tration  from  an  integral  material  balance  for  the  sulfuric 
acid.  The  discharge  process  was  simulated  and  the 
discharge  time  was  calculated  as  a  function  of  the  plate 
distance.  The  discharge  time  was  taken  to  be  completed 
when  the  porosity  reached  0.1  at  any  point.  Their  results 
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showed  that  a  certain  optimum  plate  distance  exists 
at  which  the  cell  capacity  is  highest. 

A  similar  extension  of  the  porous  electrode  models 
to  a  multi-layer  cell  was  made  by  Tiedemann  and 
Newman  [2].  Sunu  [3]  presented  a  one-dimensional 
model  for  a  lead/acid  cell  consisting  of  a  negative 
electrode,  a  porous  separator,  a  reservoir  of  stagnant 
electrolyte  and  a  positive  electrode.  The  model  was 
verified  by  comparison  of  the  predicted  and  the  ex¬ 
perimentally  measured  cell  voltages  as  a  function  of 
time  for  different  discharge  currents.  White  et  al.  [4] 
developed  a  model  similar  to  that  of  Sunu,  which  was 
later  extended  to  the  starved  and  the  hermetically  sealed 
[5,6]  lead/acid  cell.  Recently  Bernardi  et  al.  [7]  have 
treated  the  two-dimensional  case. 

In  the  above  papers,  the  models  were  mainly  used 
for  simulating  galvanostatic  discharge  of  the  cell.  Re¬ 
cently  Maja  et  al.  [8]  presented  a  model  that  was  also 
used  for  simulating  fast  charging  of  lead/acid  batteries. 

In  this  paper,  a  model  has  been  derived  for  a  lead/ 
acid  cell  consisting  only  of  a  positive  and  a  negative 
electrode  mechanically  compressed  with  a  porous  sep¬ 
arator  in  between  them.  This  cell  design  has  proved 
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to  yield  a  higher  life  under  electric-vehicle  conditions 
than  the  ordinary  starter  battery-type  [9], 

The  previous  models  mentioned  above  use  the  trans¬ 
port  equations  for  concentrated  binary  electrolytes.  The 
rate  equations  for  the  charge-transfer  processes  are 
generally  assumed  to  follow  an  ordinary  Butler-Volmer 
equation,  and  the  values  of  the  kinetic  parameters  have 
generally  not  been  critically  evaluated  or  measured. 
The  electrochemically  active  surface  area  is  assumed 
to  decrease  linearly  with  the  local  extent  of  discharge 
according  to  a  relation  that  does  not  take  into  account 
the  strong  influence  of  the  current  density  on  the  active 
mass  utilization.  In  the  present  paper  a  special  emphasis 
has  been  put  on  using  experimentally  determined  kinetic 
equations  and  relations  for  the  influence  of  structural 
changes,  as  well  as  on  the  experimental  validation  of 
the  model.  A  simple  method  for  the  numerical  solution 
of  the  coupled  differential  equations  for  the  different 
regions  will  also  be  presented. 


2.  Mathematical  model 

The  cell  investigated  is  schematically  described  in 
Fig.  1.  It  consists  of  a  porous  lead  dioxide  electrode, 
a  porous  negative  lead  electrode  and  a  porous  separator 
region  between  both  electrodes.  This  cell  design  with 
no  free  electrolyte  between  the  electrodes,  which  are 
compressed  with  separators  in  between  them,  is  of 
interest  for  obtaining  a  long  cycle  life  of  the  cell  [9]. 
In  both  the  electrodes  the  current  is  collected  by  a 
lead  grid.  In  the  cross  section  shown  in  Fig.  1  only 
two  horizontal  bars  of  the  grid  is  shown  in  the  positive 
electrode. 

The  basic  assumptions  underlying  the  model  are 
essentially  the  same  as  in  earlier  papers  (2-4).  The 
separator  is  considered  as  completely  saturated  with 
the  electrolyte  solution.  It  is  assumed  that  the  sulfuric 
acid  is  completely  dissociated  into  H+  and  HS04  ~ 
ions,  respectively.  The  model  equations  can  then  be 
derived  on  the  basis  of  the  transport  equations  for  a 
concentrated  binary  electrolyte  given  by  Newman  and 
Tiedemann  [10],  with  the  volume-average  velocity  taken 
as  the  reference  velocity.  Calculations  using,  e.g.,  Eq. 


Fig.  1.  Schematic  drawing  of  the  modelled  cell  showing  the  three 
different  regions. 


(32)  in  Ref.  10,  show  that  the  convective  term  can 
actually  be  neglected  in  the  acid  balances;  the  decrease 
in  porosity,  due  to  the  greater  molar  volume  of  lead 
sulfate  than  that  of  lead  or  lead  dioxide,  is  to  a  great 
extent  compensated  by  the  decrease  in  electrolyte  vol¬ 
ume,  due  to  the  higher  molar  volume  of  sulfuric  acid 
compared  with  water.  The  model  equations  for  the 
three  different  regions  are  thus  obtained  by  neglecting 
the  convective  term  and  applying  Eqs.  (8),  (11),  (12), 
(30)  and  (34)  in  Ref.  [10].  The  equations  can  be  written 
in  a  dimensionless  form  by  introducing  the  following 
dimensionless  variables  and  parameters: 

C=c/c0 

z=x/L 


V'^VF/RT 
i=i2/I 
T=tD0/L2 
a  =  KqRT/ILF 
f=IL/2FD0c0 
8  =  SJoL/I 

Transport  parameters  with  index  0  denote  free  elec¬ 
trolyte  conditions. 

2.1.  Positive  electrode,  0<z<Lp/L 

The  electrode  reaction  during  discharge  is: 

Pb02(s)  +  3H  +  +  HS04  ~  +2e~  — ►  PbS04(s)  +  2HzO 

The  material  balance  for  sulfuric  acid  in  the  pores 
of  the  lead  dioxide  electrode  can  be  written: 


6<c_  „  a/ 

9x  dz  \  D0  dz  J  p'  J  9z 


(1) 


The  current  balance  is: 


The  local  current  density  varies  with  the  overpotential 
according  to  the  following  Butler-Volmer  equation  for 
the  electrode  kinetics: 


j  =jo.  P[exp((2  -  «c)t}')  -  exp( -acr,')]  (3) 

where  j0  p  varies  with  the  concentration. 
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2.2.  Separator,  Lp/L<z<(Lp+h)/L 


The  material  balance  is  similar  to  that  for  the  positive 
electrode  but  with  no  source  term: 


ac  _  _a_  /£>eff3C\  dt°  ac 

0r  dz  \  D0  dz)  d C  dz 


(4) 


The  equation  for  the  current  density  through  the 
separator  is: 


—  -a  2tn)dln^dC  -  — 
dz  +  d  C  dz  flKeff 


(5) 


2.3.  Negative  electrode,  (Lp+h)/L<z<l 


The  electrode  reaction  in  the  negative  electrode  is: 
Pb(s)  +  HS04-  — ►  PbS04(s)  +  H+  +  2e" 

The  material  balance  is  analogous  to  that  for  the 
positive  electrode: 


dC 

6  0T 


_0_ 

02 


v  Do  dz) 


+  (l-2t°++(Kp,n-Kr,n)c0C)/| 


An  a. 


(6) 


The  current  balance  can  be  formulated  as: 

3  (  Keff0T?'\  3  /  Xett  n  ?  n  s  d  \tl(fc)  3C\ 

7)+7z[a^(1-2,^-^Tz) 


C=  1,  X-0,  6=  q,  for  0<z<  1 

(9) 

(ii)  for  z  =  0  and  z  =  l: 

^=0  and  5^=0 

3z  3z 

(10) 

(iii)  for  z=Lp/L: 

=  -(3-2^°  +2  —  d  ln(fc)  dC 

dz  okcH  \  +  c„  )  dC  dz 


(iv)  for  z-(Lp+H)/L: 


D-* 


=Dlft 


ac 

dz 


(11) 


and, 

^  =  _!Sl_  _/i  _2t°  1  d  ln(^c) 
dz  flKeff  ^  +)  d  C  dz 


(12) 


In  addition,  the  condition  that  the  concentration  is 
a  continuous  function  throughout  the  cell  has  been 
used.  The  boundary  conditions  in  Eqs.  (11)  and  (12) 
connect  the  concentration  profiles  in  the  three  different 
regions.  The  cell  voltage  can  be  calculated  by  means 
of  the  following  equation: 

£«n-£S.  (C(LP))+  y 


"  So,  n  jo,  n 


(7) 


The  rate  of  the  electrochemical  reaction  is  [11]: 


l-exp(27?') 

| JT -exp (act)') 


(8) 


Although  the  notations  for  transport  parameters, 
porosity  and  other  parameters  are  the  same  for  the 
different  regions,  their  values  are  different  in  the  dif¬ 
ferent  regions. 

In  the  equations  above,  the  electrolyte  potential  has 
been  defined  versus  a  hypothetical  lead  dioxide/lead 
sulfate  reference  electrode  in  the  lead  dioxide  electrode 
and  versus  a  lead/lead  sulfate  reference  electrode  in 
the  separator  and  in  the  negative  electrode.  These 
different  potentials  are  therefore  defined  only  within 
their  respective  region.  The  total  cell  voltage  can  be 
calculated  by  means  of  Eq.  (13),  see  below. 

The  following  initial  and  boundary  conditions  can 
be  formulated: 

(i)  for  f  =  0: 


X  [tj(Lp)  -  rj(Lp  T/i)  +  <)>(LP)  -  <t>(Lp +h)]  (13) 

where  E^n(C)  is  the  equilibrium  cell  voltage  at  con¬ 
centration  C.  This  value  can  be  found  in  reference 
tables.  The  approximate  polynomial  function  formulated 
by  White  and  co-workers  [5]  has  been  used  in  this 
work. 

The  effects  of  structural  changes  have  been  taken 
into  account  by  means  of  the  parameter  X,  the  degree 
of  discharge: 


Xi  = 


f  ^  dr 
D0q„  ,  J  dz 


(14) 


The  varying  porosity  in  the  porous  electrodes  can 
be  calculated  as  a  function  of  X: 

e=e0,l-ki(l-e0i)Xi  (15) 

Effective  transport  parameters  are  assumed  to  vary 
according  to  the  following  relationship,  written  for  the 
effective  conductivity  as  an  example: 

<c) 


(16) 


220 


J.  Landfors  et  al.  I  Journal  of  Power  Sources  55  (1995)  217-230 


where  A,  is  the  ratio  between  the  effective  conductivity 
in  the  active  mass  in  its  fully  charged  state  and  the 
conductivity  in  free  electrolyte.  k(C)  is  the  conductivity 
of  the  free  electrolyte  at  concentration  C.  A  similar 
relation  can  be  formulated  for  the  effective  diffusivity. 

The  active  surface  area  decreases  during  discharge 
according  to  a  relation  that  has  borrowed  its  essential 
features  from  the  empirical  Peukert  equation  [11]: 

<i7> 

This  relation  takes  into  account  the  important  fact  that 
the  electrode  surface  will  be  covered  more  rapidly  with 
smaller  lead  sulfate  crystals  when  the  current  density 
is  higher. 


3.  Influence  of  the  grid  and  its  shape 

The  presence  of  the  current-collecting  grid  is  difficult 
to  take  into  account  in  the  model.  The  grid  occupies 
about  15  to  20%  of  the  volume  of  the  electrode. 
Furthermore,  the  shape  of  the  grid  is  such  that  its 
fraction  eg  of  the  vertical  cross  section  varies  linearly 
from  zero  at  the  exterior  surface  to  a  maximum  at  the 
midplane  of  the  electrode,  see  Fig.  1.  An  exact  analysis 
would  require  a  two-dimensional  treatment.  Here  two 
different  simplified  approaches  will  be  made,  both  lead¬ 
ing  to  a  one-dimensional  analysis. 

In  the  first  case,  the  grid  cross  section  is  taken  to 
be  rectangular,  constituting  a  constant  fraction  of  the 
plate  cross  section,  so  that  the  problem  becomes  a 
truly  one-dimensional  one.  The  true  geometric  current 
density  and  effective  transport  parameters  have  to  be 
calculated  from  the  values  measured  with  respect  to 
the  total  geometric  area  including  the  grid  by  dividing 
them  by  the  factor  (l-eg),  where  eg  is  the  average 
fraction  of  the  vertical  cross  section  that  the  lead  grid 
occupies.  In  this  case  the  theoretical  charge  output  per 
unit  volume  of  active  mass  should  be  written: 


Table  1 

Values  of  the  parameters  used  in  the  calculations;  temperature  = 
23  °C 


Parameter 

Value  or  expression 

Ref. 

Co 

5.0  mol  dm-3 

D„ 

2.9 x  10-9  m2  s-1 

[4] 

DJD0 

0.706  +  0.294C 

[14] 

A, 

(VP-V,)/VTJ 

lp 

1.9  mm 

measured 

L„ 

1.7  mm 

measured 

L 

7.0  mm 

measured 

d  ln(/c)  , 

4.4763/, +  9.605/2 

[12] 

dC 

(6-/2) 

0.8224  -  0.0725C-0.0302C  2  or  different 

constant 

Soi°0 

25200  (p),  1.9X105  (n)  A  m"3 

[13,11] 

Jin, 

-1X105  A  m~3 

[11] 

Kp.  p,  VPf„ 

48.9X10-3  m3  kmol-1 

[14] 

vr,P 

25.5  X 10-3  m3  kmor' 

[14] 

K.n 

18.3 X 10-3  m3  kmol-1 

[14] 

Ap 

0.19 

[13] 

A„ 

0.33 

[11] 

A, 

varies  according  to  Figures 

n 

1.35  (p),  1.4  (n) 

<*c 

1.5  (p),  0.9  (n) 

[13,11] 

eg 

0.17  (average) 

measured 

«o 

0.55-0.58  (p),  0.60-0.66  (n) 

measured 

*0 

76  n-1 

[4] 

* 

«o(0.20 + 2. 1C  —  1.3C2)  (0.2<C<  1) 
«o2.84C  (0<C<0.2) 

[14] 

•Where/,  =0.07941  exp(2.9842C)  and /2  =  0.06645  exp(-6.4033C). 


eg=4egzf-  (20) 


(ii)  Lp/2L<z<Lp/L: 


(iii)  (Lp+/t)/L<z<(Lp+/t+Ln/2)/L: 


q‘=2F  (1  — e°.')  (18> 

In  the  second  case,  the  macrohomogeneous  approach 
is  extended  to  include  also  the  grid  material  in  such 
a  way  that  the  total  effective  porosity  in  terms  of 
electrolyte  volume  per  electrode  volume  on  the  left- 
hand  side  in  the  acid  balances  is  written: 

e=em(l-eg)  (19) 

where  em  is  the  porosity  of  the  active  material.  For 
the  two  different  electrodes  eg  varies  according  to  the 
following  relations: 

(i)  0^z<Lp/2L: 


(iv)  (Lp+/i+Ln/2)/L<z<  1: 


The  geometric  current  density  and  mass  fluxes  are 
defined  with  respect  to  the  total  cross-sectional  area, 
including  the  grid. 

The  transport  parameters  must  in  this  case  be  mod¬ 
ified  according  to  the  following  relation  for  the  effective 
diffusion  coefficient,  compare  with  Eq.  (16): 

(l-fg)D(C)  (24) 
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The  theoretical  charge  output  per  unit  electrode 
volume  and  the  initial  active  surface  area  also  vary 
along  the  depth  of  the  electrode: 

qi  =  2F^(l-e0J (l-eg)  (25) 

S0.,  =  SUl-eg)  (26) 

Input  parameters  are  given  in  Table  1.  Values  of 
the  transport  parameters  in  free  electrolyte  have  been 
taken  from  data  in  the  literature.  The  mean  activity 
coefficient  for  the  sulfuric  acid  dissociated  as  a  binary 
electrolyte,  has  been  calculated  from  a  correlation  of 
data  for  the  mean  activity  coefficient  of  hydrogen  and 
sulfate  ions  at  complete  dissociation  of  sulfuric  acid 
[12],  The  conductivity  and  diffusivity  of  5.0  M  sulfuric 
acid  were  taken  from  values  given  for  25  °C  in  Ref. 
[1],  which  were  recalculated  to  23  °C  using  the  equations 
for  temperature  dependence  for  these  parameters  given 
in  the  same  paper. 

Kinetic  parameters  for  the  lead  and  lead  dioxide 
electrodes  were  taken  from  Refs.  [11]  and  [13]  re¬ 
spectively.  In  the  latter  work  fitted  values  of  the  cathodic 
charge-transfer  coefficient  varied  around  1.5.  Using 
instead  a  value  of  2,  as  found  for  other  lead  dioxide 
electrodes  used  in  previous  works  [14,15]  did  not  sig¬ 
nificantly  change  the  results. 


4.  Numerical  procedure 

The  system  of  non-linear  differential  equations  was 
solved  by  means  of  the  finite  difference  method.  Within 
each  region  (positive  electrode,  separator  and  negative 
electrode,  respectively)  a  uniform  grid  was  chosen.  The 
first  derivative  dC/Sz  in  the  grid  point  with  number  i 
was  approximated  as  ( Ci+1-Ci_1)/(2  A h),  where  Ah 
is  the  step  length.  This  gives  a  second-order  approx¬ 
imation  of  the  first  derivatives.  The  second  derivatives 
like,  for  example,  the  first  term  on  the  right-hand  side 
of  Eq.  (1),  cannot  be  approximated  that  easily.  The 
reason  for  this  is  that  the  effective  diffusion  coefficient 
depends  also  on  concentration  and  porosity,  i.e.,  varies 
with  the  coordinate  z.  Therefore,  the  finite  difference 
approximation  at  the  point  with  number  i  can  be  taken 
in  the  form: 


(27) 

The  value  of  the  diffusion  coefficient  at  the  point 
(f  + 1/2)  is  calculated  in  terms  of  the  value  of  concen¬ 
tration  at  this  point,  Ci+in,  which  can  be  estimated  as 
(C,  +  j  +  C,)/2.  The  value  of  the  porosity  in  a  point 
between  two  mesh  points  is  calculated  in  a  similar  way. 


Eq.  (27)  gives  an  approximation  with  a  precision  of 
somewhat  less  than  second  order,  depending  on  how 
significantly  the  diffusion  coefficient  varies. 

The  first-order  accuracy  implicit  stepping  technique 
has  been  used  for  discretization  of  the  time  derivatives. 
The  approximation  of  boundary  conditions  is  very  im¬ 
portant.  It  is  well  known  that  the  boundary  condition 
dCdz= 0,  written  in  the  simplest  form: 

Ca-Co=0  (28) 

gives  only  a  first-order  accuracy,  which  would  make 
the  second-order  approximation  for  the  differential 
equation  lose  its  attraction.  White  [5,6,16]  proposed  to 
use  a  three-point  approximation  for  each  of  the  first 
derivatives  at  the  boundary  between  two  regions.  This 
gives  an  equation  for  the  boundary  condition  that 
involves  five  mesh  points.  An  algorithm  for  solving  such 
a  problem  was  proposed  in  Ref.  [17]. 

In  this  work  we  used  instead  the  approximation  of 
the  boundary  condition  which  is  standard  in  fluid  dy¬ 
namic  numerical  works.  The  mesh  grid  was  shifted  so 
that  all  boundaries  were  half  way  between  two  grid 
points,  and  imaginary  points  were  added  at  a  half  step 
outside  the  regions.  In  this  case,  Eq.  (28),  where  now 
i=0  corresponds  to  an  imaginary  point,  has  second- 
order  accuracy. 

For  the  boundary,  the  value  of  the  effective  diffusion 
coefficient  was  taken  as  the  harmonic  mean  value  of 
the  values  of  the  effective  diffusion  coefficient  in  the 
neighbouring  grid  points  in  the  electrode(e)  and  sep¬ 
arators),  respectively: 


The  same  procedure  can  be  used  for  the  equations 
for  the  electric  potential,  Eqs.  (2)  and  (7).  No  interface 
conditions  for  the  concentration  are  required.  Three 
equations  for  the  concentration  in  the  three  different 
regions,  Eqs.  (1),  (4)  and  (6),  are,  in  fact,  made  into 
one  quasi-linear  equation  with  variable  coefficients.  On 
the  boundaries  of  the  separator  these  coefficients  change 
stepwise,  but  from  physical  and  numerical  points  of 
view  it  remains  the  same  equation.  In  this  way  we  get 
a  standard  tridiagonal  matrix  for  the  full  system,  in¬ 
dependently  of  the  number  of  the  inner  boundaries. 

The  algorithm  used  for  solving  the  discretized  equa¬ 
tions  is  described  below: 

(i)  Eq.  (2)  for  the  potential  within  the  positive 
electrode  with  boundary  conditions  given  in  Eqs.  (10a) 
and  (lib)  was  linearized  and  solved  by  an  over-relaxation 
procedure.  For  the  solution  of  implicit  space  derivatives 
a  simple  Thomson  algorithm  was  used.  In  order  to 
have  a  stable  convergence  from  arbitrary  ‘initial’  values, 
an  artificial  viscosity  term  was  added.  At  the  end  of 
this  stage  two  conditions  were  checked:  (i)  a  stable 
potential  profile  has  been  obtained,  and  (ii)  the  boundary 
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conditions  given  in  Eqs.  (10a)  and  (lib)  for  the  full 
equations  (not  only  for  the  linearized)  have  been  fulfilled 
with  sufficient  accuracy  (about  10-8). 

(ii)  The  changes  in  porosity  and  active  surface  at 
this  potential  distribution  is  calculated  from  Eqs.  (14) 
to  (17).  The  new  values  of  porosity  requires  a  correction 
of  the  potential  profile,  i.e.,  returning  to  the  step  (1). 
Then  we  can  correct  porosity  one  more  time.  The 
procedure  according  to  the  steps  (1)  to  (2)  converged 
very  fast,  it  took  only  2  to  3  iterations. 

(iii) ,  (iv)  An  analogous  procedure  for  the  potential 
and  porosity  distributions  within  the  negative  electrode 
must  be  done.  Note,  that  we  do  not  solve  the  equation 
for  the  potential  in  the  separator.  Eq.  (5)  for  the 
potential  distribution  can  be  integrated  with  respect  to 
the  coordinate.  Therefore,  if  we  know  the  concentration 
profile  within  the  separator,  we  can  find  the  potential 
distribution  immediately. 

(v)  The  overvoltage  and  porosity  profiles,  obtained 
in  steps  (i)  to  (iv)  were  substituted  into  the  equation 
for  concentration.  The  single  equation  for  the  complete 
multi-region  system  was  used.  This  equation  has  also 
non-linear  factors,  depending  upon  the  concentration. 
Therefore,  an  iteration  procedure  is  necessary.  After 
calculations  of  new  values  of  concentration,  the  new 
potential  distribution  was  calculated  from  steps  (i)  and 
(iii)  and  then  porosity  corrections  were  made,  steps 
(ii)  and  (iv). 

The  main  convergence  loop  including  steps  (i)  to  (v) 
gives  in  general  good  results  already  after  two  or  three 
iterations,  but  in  some  complicated  cases,  for  example, 
at  the  end  of  discharge,  when  a  region  with  low  con¬ 
centration  appears,  the  number  of  iterations  can  grow 
significantly. 

At  the  end  of  procedure  (i)  to  (v),  we  obtained  the 
solution  of  the  non-linearized  system  on  the  next  time 
step,  using  implicit  time-derivative  technique.  After  that, 
the  next  time  step  can  be  taken. 

The  main  advantages  of  this  algorithm  are: 

(i)  It  is  absolutely  stable,  it  converged  for  any  time 
step  and  length  coordinate  step,  as  well  as  for  small 
values  of  the  concentration.  In  the  case  when  conver¬ 
gence  becomes  bad,  the  program  automatically  reduces 
the  time  step.  It  gives  a  possibility  to  calculate  up  to 
very  high  degrees  of  discharge  when  other  algorithms 
tried  could  not  work. 

(ii)  It  gives  us  the  possibility  to  check  the  accuracy 
of  the  calculations  of  each  component  separately  and 
of  the  system  as  a  whole. 

(iii)  At  the  same  time,  we  do  not  loose  the  attractive 
feature  of  the  Thomson  algorithm,  namely  linearity; 
the  calculation  time  is  proportional  to  the  number  of 
grid  points  (rather  than  to  the  square  of  this  number). 

(iv)  A  discretization  of  a  second-order  accuracy  was 
used  for  the  equations  as  well  as  the  interface  and 
boundary  conditions. 


The  execution  time  on  a  personal  computer  is  several 
minutes,  which  allows  us  to  observe  the  system  evolution 
directly  on  the  screen. 

Further  details  about  the  stability  and  conservative 
properties  of  the  numerical  method  are  given  in  the 
Appendix. 


S.  Experimental 

The  electrodes  used  in  this  study  were  pasted  elec¬ 
trodes  for  commercially  available  SLI-batteries  man¬ 
ufactured  by  Tudor  AB  at  Nol,  Sweden.  The  full  size 
dimensions  were  115  mm  X 155  mm  with  a  thicknesses 
of  approximately  1.7  and  1.9  mm  for  the  lead  and  the 
lead  dioxide  electrodes,  respectively.  The  current-col¬ 
lector  grids  were  made  of  a  lead  alloy  with  2.5%  Sb. 
The  experiments  were  performed  on  electrode  samples 
that  had  been  cut  from  the  full  size  plates  to  about 
51  mm  in  height  and  34  mm  in  width.  The  samples 
were  taken  both  from  the  upper  and  the  lower  parts 
of  the  plates,  so  that  one  side  of  each  electrode  sample 
consisted  of  the  grid  bar  to  which  a  3  mm  lead  rod 
was  welded  as  the  current  collector. 

The  electrolyte  used  during  the  preparative  cycling 
and  before  discharge  experiments  was  in  all  cases  5.05 
M  sulfuric  acid  prepared  from  98%  sulfuric  acid  (Merck, 
pro  analysi)  and  purified  water  from  a  combined  Milli- 
RO  15  and  Milli-Q  water  purification  system  from 
Millipore.  The  concentrations  were  determined  by 
acid-base  titrations  and/or  refractive  index  measure¬ 
ments. 

Lead  and  lead  dioxide  electrodes  were  precycled  in 
pairs  in  glass  beakers;  the  electrodes  were  separated 
by  a  microporous  separator  to  avoid  short-circuiting. 
The  electrodes  were  allowed  to  soak  in  the  electrolyte 
for  about  1  h  before  they  were  automatically  cycled  at 
a  constant-current  density  of  100  Am-2  according  to 
the  following  scheme: 

1.  charge  to  oxygen  evolution; 

2.  rest  for  3600  s; 

3.  discharge  for  14605  s,  corresponding  to  approximately 
80%  depth-of-discharge  at  rated  capacity  of  the  cell; 

4.  rest  for  3600  s; 

5.  recharge  for  17  526  s,  corresponding  to  20%  over¬ 
charge,  and 

6.  steps  2  to  5  were  repeated  in  total  three  times. 
After  the  preparative  cycling,  the  electrolyte,  in  the 

beakers,  was  changed  to  ensure  a  concentration  of  5.05 
M.  The  beakers,  with  electrolyte  and  electrodes,  were 
placed  under  vacuum  maintained  by  a  water  jet  for  a 
period  of  at  least  10  h.  In  order  to  remove  as  much 
as  possible  of  the  gas  in  the  porous  electrodes  and  to 
replace  the  gas  with  the  electrolyte.  After  evacuation, 
the  electrolyte  was  changed  once  again.  The  temperature 
during  the  discharge  experiments  was  kept  at  23  °C. 
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The  so-prepared  electrodes  were  mounted  into  an 
experimental  cell  of  which  a  schematic  drawing  is  given 
in  Fig.  2.  The  cell  consisted  of  two  thick  endplates 
made  of  poly(tetrafluoroethylene)  (PTFE)  with  a  U- 
formed  distance  made  of  poly(methyl  methacrylate) 
(PMMA)  forming  a  well-defined  cell.  The  total  distance 
over  the  electrodes  was  7.0  mm.  The  height  and  the 
width  of  the  cell  matched  the  form  and  dimensions  of 
the  electrodes.  In  each  endblock  there  was  a  hole 
perpendicular  to  the  electrode  surface.  The  hole  was 
used  to  measure,  via  a  capillary,  the  electrode  potential 
at  the  rear  of  each  electrode  by  Hg/Hg2S04  reference 
electrodes  in  saturated  K2S04  manufactured  by  Ra¬ 
diometer.  Three  glass-fibre  separators  (BG  200  17  by 
Hollingsworth  and  Vose),  with  the  same  dimensions 
as  the  electrodes,  were  placed  between  the  electrodes. 
The  separators  were  immersed  in  sulfuric  acid  and 
evacuated  in  the  same  way  as  the  electrodes.  A  me¬ 
chanical  pressure  was  applied  perpendicularly  to  the 
electrode  plates  to  seal  the  cell  and  to  get  a  well- 
defined  electrode  distance.  Any  free  acid,  i.e.,  acid  that 
was  not  soaked  by  the  electrodes  or  the  separator,  was 
removed  from  the  cell  container  with  a  syringe.  This 
setup  defined  the  cell  at  time  equal  to  zero. 

The  cell  was  then  galvanostatically  discharged  at 
different  current  densities  varying  from  100  to  1000  A 
m  ~ 2.  During  discharge,  the  cell  voltage  and  the  electrode 
potentials  at  the  rear  of  the  electrodes  were  measured 
on  x-t  recorders.  The  discharge  was  stopped  at  a  cell 
voltage  of  1.0  V.  The  cell  was  then  immediately  di¬ 
sassembled.  The  remaining  amount  of  sulfuric  acid  was 
then  separately  leached  out  from  each  cell  component 
by  soaking  in  multiple  small  portions  of  purified  water 
until  the  leaching  water  was  neutral.  The  amount  of 
acid  in  each  cell  component  could  then  be  determined 
by  an  acid-base  titration.  Some  cells  were  disassembled 
without  discharge  in  order  to  determine  the  amount 
of  acid  and  thereby  the  effective  porosity  of  the  elec¬ 
trodes  and  the  separators.  The  concentration  profiles 
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Fig.  2.  Schematic  drawing  of  the  cell  used  in  the  discharge  experiments. 


in  the  cells  after  discharge  could  thus  be  calculated 
from  the  remaining  amount  of  acid  according  to  the 
titrations  as  well  as  the  final  porosity  which,  for  the 
two  electrodes,  could  be  calculated  from  the  initial 
porosity  and  the  change  in  porosity  due  to  the  formation 
of  lead  sulfate  from  lead  and  lead  dioxide,  respectively. 

The  distribution  of  lead  sulfate  along  the  thickness 
of  the  electrodes  was  determined  with  an  electron- 
probe  microanalyser  (ARL  SEMQ  42)  by  recording  the 
S  Ka  radiation.  The  electrodes  for  the  sulfate  analyses 
were  leached  in  acetone  instead  of  water.  The  electrodes 
were  then  dried  at  50  °C  and  representative  samples 
from  the  central  part  of  the  electrodes  were  taken  out. 
The  samples  were  then  moulded  into  an  epoxy  resin 
under  vacuum,  polished  with  diamond  paste  and  sput¬ 
tered  with  a  conductive  layer  of  carbon  before  analysis. 


6.  Results  and  discussion 

6.1.  Effects  of  different  parameters 

The  computer  program  was  first  tested  with  the  model 
equations  and  input  parameters  used  by  Nguyen  and 
White  in  Ref.  [16].  The  results  obtained  were  in  good 
agreement  with  those  of  Fig.  2  in  that  work.  Similar 
results  obtained  with  the  model  and  input  parameters 
used  in  this  work  are  shown  in  Fig.  3  for  a  current 
density  of  1000  Am-2.  This  Figure  shows  that  the 
discharge  capacity  at  a  high  rate  of  discharge  is  mainly 
determined  by  an  acid  depletion  in  the  positive  elec¬ 
trode. 


Distance 


Fig.  3.  Model  predictions  of  the  changes  in  sulfuric  acid  concentration 
in  a  cell  at  a  discharge  with  1000  A  m~2.  Constant  /+  =0.79.  Porosity 
of  positive  plate = 0.58,  of  negative  plate  =  0.66,  and  of  separator  =  0.53. 
Thickness  of  positive  plate  =  2.15  mm,  of  negative  plate  =  1.8  mm, 
and  of  separator  region  =  2.4  mm.  ATj  =  15  X  107  (A  m~2)135  s  for 
both  electrodes. 
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In  these  calculations  a  constant  value  of  the  trans¬ 
ference  number  was  used.  It  turned  out  that  calculations 
with  a  varying  transference  number  did  not  give  con¬ 
centration  profiles  that  were  in  agreement  with  the 
acid  conservativity  condition.  The  calculated  consump¬ 
tion  of  sulfuric  acid  could  deviate  by  up  to  10%  from 
the  value  calculated  by  means  of  an  overall  acid  balance 
using  Faraday’s  law.  All  subsequent  calculations  were 
therefore  performed  with  a  constant  transference  num¬ 
ber. 

Initial  calculations  showed  that  the  shape  of  the  grid 
that  was  assumed  in  the  model  did  not  affect  the  results 
significantly.  The  main  difference  was  that  the  porosity 
profiles  calculated  by  means  of  using  Eqs.  (20)  to  (23) 
exhibited  a  slight  maximum  in  the  interior  of  the 
electrode,  see  Fig.  5.  The  calculations  reported  here 
were  performed  using  Eqs.  (20)  to  (23). 

The  effects  of  the  parameters  K'a  and  n  are  shown 
in  Figs.  4  and  5.  When  n  =  1  and  the  maximum  possible 
utilization  of  the  active  electrode  material  is  equal  to 
the  theoretical  value,  then  mainly  the  outer  layer  of 
the  electrode  has  been  discharged  and  the  local  degree 
of  discharge  has  its  maximum  at  the  outer  surface  of 
the  electrode.  With  n  =  1.35  and  a  value  of  K'a,  which 
corresponds  to  a  maximum  degree  of  utilization  that 
is  lower  than  the  theoretical  one,  the  maximum  of  the 
current-density  distribution  moves  progressively  inwards 
as  the  outer  layers  of  the  electrode  are  de-activated 
by  the  formed  lead  sulfate.  At  a  constant  value  of  n 
a  larger  value  of  K'a  means  in  general  a  higher  possible 
utilization  of  the  active  electrode  material  and,  there- 


Fig.  4.  Predicted  concentration  profiles  for  different  values  of  K'a 
for  the  positive  electrode:  (dotted  line)  K'„  =  3X 107  (A  m-2)1’35  s; 
(bold,  solid  line)  X’j  =  7xl07  (A  m-2)1-35  s  and  (dashed  line) 
Kj  =  20xl07  (A  m-2)1-35  s.  Other  data  as  in  Fig.  3.  Dot-dashed  line 
for  n  =  l  instead  of  1.35  and  only  theoretical  limitation  of  the  local 
degree  of  discharge.  Discharge  time  =  472,  700,  750  and  760  s, 
respectively. 


Distance 

Fig.  5.  Model  predictions  of  the  changes  in  porosity  for  the  different 
cases  in  Fig.  4. 


fore,  a  higher  discharge  capacity.  A  lower  value  of 
K'a  for  the  positive  electrode  forces  the  reaction  deeper 
into  the  interior  of  the  electrode  due  to  a  more  rapid 
de-activation  of  the  parts  of  the  electrode  that  are  close 
to  the  interface.  For  n>l,  the  high  initial  current 
density  in  the  outermost  layer  of  the  electrode  causes 
the  electrode  material  in  this  layer  to  become  de¬ 
activated  at  a  relatively  low  degree  of  discharge.  The 
degree  of  discharge  therefore  has  a  maximum  at  some 
distance  from  the  interface  between  the  electrode  and 
the  separator.  This  is  clearly  pronounced  for  K'a  =  7  x  107 
(A  m-2)1-35  s  for  the  positive  electrode  in  Fig.  5.  This 
shape  of  the  curve  for  the  degree  of  discharge  is  in 
better  agreement  with  the  distributions  of  formed  lead 
sulfate  measured  with  electron-probe  microanalysis.  It 
was  found  that  a  reasonable  value  of  K'a  lies  in  the 
region  of  (7-15) X107  (A  m-2)135  s  for  the  positive 
electrode.  In  this  region  the  discharge  is  mainly  limited 
by  an  acid  depletion  in  the  positive  electrode. 

In  contrast,  for  the  lowest  value  of  K'a  used  in  Figs. 
4  and  5,  the  maximum  local  mass  utilization  is  very 
low  and  the  discharge  capacity  becomes  unrealistically 
low.  Figs.  4  and  5  show  that  in  this  case  the  discharge 
capacity  would  be  limited  by  the  low  (and  fairly  uniform) 
active  mass  utilization. 

Fig.  6  shows  the  effect  of  the  separator  porosity  on 
the  discharge  capacity  at  a  high  rate  discharge  (1000 
Am-2).  It  can  be  seen  that  this  parameter  has  a  rather 
weak  influence  on  the  discharge  capacity,  although  at 
given  dimensions  of  the  different  regions  a  higher 
separator  porosity  means  both  a  higher  effective  dif- 
fusivity  and  a  larger  amount  of  the  stored  acid.  The 
explanation  for  this  weak  dependence  at  high  rates  of 
discharge  is  that  mainly  the  acid  in  the  pores  of  the 
electrode  is  being  consumed  in  the  electrode  reactions 
during  the  relatively  short  discharge  process,  which  is 
terminated  by  the  depletion  of  sulfuric  acid  in  a  zone 
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Fig.  6.  Predicted  concentration  profiles  for  different  values  of  the 
separator  porosity:  (solid  line)  0.96,  (dashed  line)  0.7  and  (lower 
line)  0.5.  A,  for  the  separator  was  taken  as  eIJ.  Constant  I°+  =0.79. 
ATi  =  15xl07  for  both  electrodes.  Discharge  time  =  724,  719  and  699 
s,  respectively. 

just  behind  the  de-activated  outer  layer  of  the  positive 
electrode. 

For  low  rates  of  discharge,  when  the  discharge  capacity 
is  higher  and  the  consumption  rate  of  sulfuric  acid  in 
the  electrodes  can  more  easily  be  balanced  by  the 
transport  rate  of  acid  from  the  separator  into  the 
electrodes,  a  high  separator  porosity  has  a  much  more 
positive  influence. 

A  similar  dependence  holds  for  the  discharge  capacity 
as  a  function  of  the  thickness  of  the  separator. 
Figs.  7  and  8  show  discharge  curves  for  a  cell  with  1.0 
mm  thick  electrodes  and  with  separators  of  different 
thicknesses.  At  a  high  rate  of  discharge  (Fig.  7)  the 
discharge  capacity  is  fairly  independent  of  the  separator 
thickness.  Surprisingly,  the  discharge  capacity  is  even 
larger  for  the  thinner  separator.  This  can,  at  least  partly, 
be  explained  by  the  higher  ohmic  potential  drop  in  the 
thicker  separator. 

For  a  lower  rate  of  discharge  (Fig.  8)  the  discharge 
capacity  increases  with  the  separator  thickness,  since 
the  acid  in  the  pores  of  the  separator  is  consumed  to 
a  higher  extent  in  this  case.  At  the  beginning  of  discharge, 
the  cell  voltage  is  the  highest  for  the  thinnest  separator, 
due  to  a  lower  ohmic  potential  drop  across  the  separator, 
while  the  relation  is  reversed  with  time  when  the  overall 
depletion  of  sulfuric  acid  becomes  determining. 

In  order  to  have  a  sufficiently  large  amount  of  elec¬ 
trolyte  in  the  case  low  rate  discharges,  the  separator 
region  must  therefore  be  made  sufficiently  thick  and/ 
or  must  have  a  sufficiently  high  porosity. 

6.2.  Comparison  with  experimental  results 

An  important  aim  of  this  study  was  to  validate  the 
model  by  comparing  results  predicted  by  the  model 


Fig.  7.  Predicted  voltage-time  curves  for  a  cell  with  1.0  mm  thick 
electrode  and  different  thicknesses  of  the  separator:  (upper  line)  1 
mm,  (line  in  between)  2  mm,  and  (lower  line)  3  mm.  Current 
density=1000  A  m-2.  10X 107  (A  m-2)135  s  for  the  positive 

electrode  and  /C,=  15X 107  (A  m-2)14  for  the  negative  electrode. 
Porosity  of  separator  =  0.804.  A,  for  the  separator  =  0.33.  Porosity  of 
negative  electrode  =  0.66.  t\  =0.79. 


Fig.  8.  As  in  Fig.  7  but  at  a  discharge  current  density  of  100  A 

with  those  obtained  experimentally.  Although  the  elec¬ 
trodes  did  not  differ  very  much  with  respect  to  measured 
porosity  and  geometric  dimensions,  and  the  preparative 
cycling  was  the  same,  the  discharge  capacity  varied 
unexpectedly  much.  This  is  shown  in  Fig.  9,  in  which 
calculated  and  measured  concentration  profiles  at  a 
high  current  density,  490  Am-2,  are  compared.  The 
experimental  discharge  capacity  in  four  different  ex¬ 
periments  varied  between  2121  and  2548  s.  The  cal¬ 
culated  discharge  time  of  2012  s  is  in  good  agreement 
with  the  lowest  value  obtained  experimentally,  which, 
on  the  other  hand,  is  8%  lower  than  the  average  and 
17%  lower  than  the  highest  discharge  time  obtained 
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Fig.  9.  Comparison  of  the  (solid  lines)  predicted  and  (dashed  lines) 
measured  concentration  profiles  at  a  discharge  with  490  A  m-2. 
r°  =0.81.  Assumed  porosity  of  negative  electrode  =  0.60,  of  positive 
electrode  =  0.58.  Porosity  of  separator  =  0.96.  A,  for  the 
separator  =  0.90.  Predicted  discharge  time  =  2012  s.  Experimental 
discharge  time  =  2121  (dashed  lines),  and  2548  s  (dotted  lines). 

experimentally.  The  measured  concentration  profiles 
are  in  fair  agreement  with  those  predicted  by  the  model. 
Both  the  experiments  and  the  theory  show  that  the 
discharge  capacity  at  high  discharge  rates  is  limited  by 
an  acid  depletion  in  the  positive  electrode. 

Similar  results  were  obtained  for  a  current  density 
of  1000  Am-2.  According  to  the  model  predictions 
in  Figs.  3  and  4,  the  discharge  time  would  be  around 
700  to  750  s,  whereas  the  experimental  results  give 
discharge  times  varying  between  830  and  940  s. 

For  200  Am-2  (Fig.  10),  the  model  predicts,  with 
input  parameters  valid  for  room  temperature,  a  dis¬ 
charge  time  that  is  4  to  7%  lower  than  that  obtained 
experimentally.  The  predicted  discharge  time  for  100 
A  m“2  (17  600  s)  (Fig.  11)  is  in  good  agreement  with 
the  experimentally  measured  discharge  time  of 
16  620-17  970  s.  The  measured  concentrations  in  the 
separator  are  somewhat  higher  than  the  predicted  val¬ 
ues,  whereas,  on  the  other  hand,  the  measured  con¬ 
centration  in  the  negative  electrode  is  again  significantly 
lower  than  the  calculated  value. 

Figs.  9  to  11  show,  that  in  comparison  with  the 
experimental  results,  the  model  generally  predicts  a 
somewhat  too  rapid  depletion  of  sulfuric  acid,  especially 
at  high  rates  of  discharge,  and  a  too  low  decrease  in 
sulfuric  acid  concentration  in  the  negative  electrode. 
A  common  cause  for  these  two  discrepancies  should 
first  be  searched  for.  One  parameter  that  affects  the 
resulting  distribution  of  acid  between  the  different 
regions  is  the  transference  number.  According  to  Eqs. 
(1)  and  (6)  a  higher  value  of  the  transference  number 
gives  a  slower  consumption  rate  of  sulfuric  acid  in  the 
positive  electrode  and  a  higher  consumption  rate  in 
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Fig.  10.  Comparison  of  the  (solid  lines)  predicted  and  (dashed  lines) 
measured  concentration  profiles  at  a  discharge  with  200  Am-2. 
Other  input  data  as  in  Fig.  9.  Predicted  discharge  time  =  6600  s. 
Experimental  discharge  time  =  6942  s. 


Fig.  11.  Comparison  of  (solid  lines)  predicted  and  (dashed  lines) 
measured  concentration  profiles  at  a  discharge  with  100  A  m-2. 
Other  input  data  as  in  Fig.  9.  Predicted  discharge  time  =17600  s. 
Experimental  discharge  time  =  17310  s. 

the  negative.  It  should  be  kept  in  mind  that  the  value 
of  this  parameter  for  the  volume-average  velocity  as 
the  reference  velocity  (as  in  this  case)  has  somewhat 
arbitrarily  been  set  equal  to  the  transference  number 
with  respect  to  the  solvent  velocity  [10],  According  to 
the  derivation  of  the  transport  equation  with  vol¬ 
ume-average  velocity  as  reference  velocity: 

f+  =  (l-c+K+)f°++c_K_f°_  (30) 

where  c+V+  and  c_K_  cannot  be  determined  indi¬ 
vidually.  Newman  has,  therefore,  arbitrarily  put 
t°+  v+  V+  =t°_  v_  V which  gives  t+=t°+.  This  relation 
has  been  used  in  our  model.  Now,  assuming  instead 
that  the  molar  volume  of  H+  =0,  which  seems  not  too 
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unrealistic,  considering  the  size  of  a  proton.  This  also 
means  that  c_V_~cV„  and  we  obtain: 

t+  =  {l-cvy++cvc  (31) 

Using  known  data  it  can  be  shown  that  t+  estimated 
from  Eq.  (31)  is  only  1%  higher  than  t°+  for  a  1  M 
sulfuric  acid  solution,  whereas  it  is  8%  higher  for  a  5 
M  sulfuric  acid  solution. 

Fig.  12  shows  the  effect  of  the  transference  number 
on  the  resulting  concentration  profile.  As  anticipated, 
the  obtained  concentration  profile  in  the  cell  is  in  better 
agreement  with  the  experimental  results  when  the  value 
of  the  transference  number  is  increased.  However,  an 
unreasonably  high  value  would  be  needed  to  give  the 
low  concentrations  in  the  negative  electrode  measured 
experimentally. 

There  are,  therefore,  probably  different  causes  for 
the  discrepancies  between  model  predictions  and  ex¬ 
perimental  results  for  the  two  different  electrodes.  In 
order  to  simulate  correctly  the  behaviour  of  the  positive 
electrode  at  high  rates  of  discharge  significantly  higher 
values  of  the  transport  parameters  are  required.  These 
higher  values  could  be  a  result  of  an  increased  tem¬ 
perature  in  the  electrode  due  to  irreversible  losses  at 
the  relatively  high  overpotentials.  This  would  also  ex¬ 
plain  the  increasing  discrepancy  between  theoretical 
and  experimental  results  with  increasing  current  density, 
since  the  thermal  effects  are  larger  at  higher  current 
densities.  If  the  thermal  effects  cause  the  unexpectedly 
high  discharge  times  at  higher  current  densities,  then 
it  is  also  reasonable  that  the  agreement  should  be  better 
at  lower  discharge  rates  when  these  thermal  effects  are 
smaller.  A  rough  estimation  of  the  adiabatic  temperature 
rise  in  the  cell  using  data  from  Ref.  [18]  gives  an 


Fig.  12.  Predicted  concentration  profiles  at  different  transference 
numbers  for  the  hydrogen  ion:  (upper  curve)  0.72,  (middle  curve) 
0.79  and  (lower  curve)  0.81.  Predicted  discharge  time  =  638,  700  and 
715  s,  respectively.  Other  data  as  in  Fig.  3  but  with  ^  =  10xl07 
(A  m-2)135  s  for  the  positive  electrode. 


average  value  of  3-5  °C  after  discharge  with  500-1000 
Am-2.  These  values  are  not  sufficiently  high  to  explain 
the  discrepancies  for  the  positive  electrode,  but  due 
to  the  highly  non-uniform  current  and  potential  dis¬ 
tributions  in  this  electrode,  the  temperature  may  rise 
to  much  higher  values  locally  in  the  critical  interface 
region  between  the  separator  and  the  positive  electrode. 
Already,  the  Joulean  heat  production  rate  in  the  outer 
layer  of  the  positive  electrode,  in  which  the  effective 
conductivity  of  the  electrolyte  is  low  due  to  a  low 
porosity  and  a  decreasing  acid  concentration,  can  be 
estimated  to  be  considerably  higher  than  the  average. 
In  order  to  investigate  the  thermal  effects  further,  a 
non-isothermal  model  incorporating  energy  balances  is 
required. 

Temperature  effects  do  not  explain  the  opposite 
deviation  in  the  negative  electrode,  in  which  the  mea¬ 
sured  sulfuric  acid  concentration  is  significantly  lower 
than  the  calculated  value.  It  should  be  stressed  that 
the  calculations  have  been  performed  with  the  lowest 
measured  value  of  the  porosity  of  the  negative  plates. 
It  is,  therefore,  not  very  likely  that  the  real  porosity 
is  actually  somewhat  lower,  which  would  otherwise 
explain  the  deviations.  The  microprobe  analyses  of  the 
sulfur  distribution  in  the  negative  electrode  after  dis¬ 
charge  show  that  the  overall  utilization  of  the  negative 
active  material  is  fairly  uniform,  whereas  the  model 
predicts  a  non-uniform  distribution  with  a  higher  degree 
of  discharge  in  those  parts  closer  to  the  electrolyte 
matrix.  This  may,  at  least  partially,  explain  why  the 
model  predicts  an  average  concentration  in  the  negative 
that  is  higher  than  the  measured  values. 


7.  Conclusions 

1.  Both  the  theoretical  and  the  experimental  results 
show  that  the  discharge  capacity  is  limited  by  the  acid 
depletion  in  the  positive  electrode  for  current  densities 
between  100  and  1000  Am'2. 

2.  At  high  rates  of  discharge,  the  amount  of  sulfuric 
acid  in  the  pores  of  the  separator  does  not  affect  the 
discharge  capacity  to  any  greater  extent.  At  lower 
discharge  rates,  on  the  other  hand,  the  discharge  capacity 
increases  strongly  with  an  increasing  amount  of  elec¬ 
trolyte  in  the  separator. 

3.  The  model  can  be  used  for  optimizing  the  geometric 
dimensions  of  the  positive  and  negative  electrodes  and 
the  separator  in  the  lead/acid  battery  for  different 
applications. 

4.  Experimental  and  theoretical  results  are  in  good 
agreement  at  medium  to  low  current  densities,  whereas 
the  model  predicts  a  somewhat  too  low  discharge  ca¬ 
pacity  at  high  current  densities,  possibly  because  it  does 
not  take  into  account  the  thermal  effects. 
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5.  Further  investigations  are  required  to  explain  the 
discrepancies  between  the  model  predictions  and  the 
experimental  results  for  the  decrease  in  the  sulfuric 
acid  concentration  in  the  negative  electrode  during 
discharge. 


8.  List  of  symbols 

a  KoRT/ILF 

c  concentration  of  sulfuric  acid  (mol  dm-3) 

C  dimensionless  concentration  of  sulfuric  acid 

D  diffusion  coefficient  of  sulfuric  acid  (m2  s-1) 
£cell  cell  voltage  (V) 

/  IL/2FD0c0 

f  activity  coefficient  in  Eqs.  (2),  (7),  (11)  and  (12) 

F  Faraday  constant  (96  487  A  s  mol-1) 

8  SofoL/I 

h  thickness  of  separator  region  (m) 

i  dimensionless  current  density  in  the  pore  elec¬ 
trolyte 

i2  current  density  in  the  pore  electrolyte  (Am-2) 

1  geometric  current  density  (A  m-2) 

j  local  true  current  density  (A  m-2) 

jfim  cathodic  limiting  current  density  (A  m-2) 
j0  exchange-current  density  (A  m-2) 

k,  (Kp>i-Fr>i)/Fr,(. 

K'zj  constant  defined  by  Eq.  (17) 

L  thickness  (m) 

M  molecular  weight  (kg  mol-1) 

n  exponent  defined  by  Eq.  (17) 

q0  i  charge  corresponding  to  the  initial  amount  of 
i  (A  s  m-3) 

R  universal  gas  constant  (8.3143  J  mol-1  K-1) 

S  specific  active  surface  area  (m-1) 

t  discharge  time  (s) 

t°+  transference  number  of  hydrogen  ion  for  solvent 
velocity  as  reference  velocity 
tT  transference  number  of  hydrogen  ion  for  vol¬ 

ume-average  velocity  as  reference  velocity 
T  absolute  temperature  (K) 

V+  molar  volume  of  the  hydrogen  ion  (m3  mol-1) 

V _  molar  volume  of  the  hydrogen  sulfate  ion  (m3 

mol-1) 

Vp  i  molar  volume  of  product  (m3  kmol-1) 

Vr  i  molar  volume  of  reactant  (m3  kmol-1) 

x  distance  ( m ) 

X  degree  of  discharge 

z  dimensionless  distance 

Greek  letters 

ac  cathodic  transfer  coefficient 

e  porosity 

eg  fraction  of  electrode  cross-sectional  area 

em  porosity  of  the  active  material 


ower  Sources 

t)  overvoltage  (V) 

17'  dimensionless  overvoltage 

k  conductivity  (fl-1  m-1) 

A  ratio  between  effective  and  free  transport  pa¬ 

rameters 

p  density  (kg  m-3) 

<f>  potential  (V) 

r  dimensionless  time 

Subscript 

0  initial  value 

e  electrode 

eff  effective  value 

n  negative  electrode 

p  positive  electrode 

s  separator 

W  water 

Superscript 

0  initial  value 
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Appendix 


Stability  and  conservative  properties  of  the 
numerical  solution 

The  algorithm  developed  by  Newman  [19]  provides 
a  stable  solution  for  a  linear  system,  independently  of 
the  variation  of  the  time  step.  For  non-linear  systems 
stability  can  be  achieved  only  for  a  finite  time  step, 
which,  as  a  rule,  can  only  be  determined  after  numerical 
experiments.  Nonetheless,  not  even  a  very  small  time 
step  can  guarantee  the  stability  of  numerical  calcula¬ 
tions.  There  are  several  causes  for  an  instability  of  the 
numerical  algorithm,  e.g.: 

(i)  Instability  due  to  the  linearization  of  non-linear 
coefficients  and  terms.  In  the  case  that  we  solve  non¬ 
linear  equations,  we  have  to  decide  up  to  what  extent 
we  are  going  to  approximate  the  non-linear  terms.  For 
example,  in  the  calculation  of  a  diffusion  coefficient 
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on  the  next  time  level  according  to  Eq.  (24),  we  have 
to  substitute  some  value  for  the  concentration  in  the 
expression  for  D(C).  It  is  evident  to  use  the  value  of 
the  concentration  from  the  previous  time  step.  Oth¬ 
erwise,  a  system  of  linear  algebraic  equations  cannot 
be  obtained.  This  gives  a  satisfactory  result  when  coef¬ 
ficients  are  changing  slowly  with  time  and  distance, 
otherwise  an  instability  will  appear.  For  example,  it 
occurs  when  we  have  a  significant  concentration  gra¬ 
dient. 

(ii)  Instability  due  to  initial  conditions.  In  our  case 
there  are  no  problems  with  initial  values  of  porosity 
and  concentration  profiles.  On  the  contrary,  the  initial 
potential  distribution  is  not  known  and  must  be  found 
from  the  solution  of  the  equation  for  the  potential 
before  the  real  calculation  can  start.  In  the  opposite 
case,  when  an  initial  potential  distribution  is  not  well 
defined,  the  values  of  the  overpotential  from  the  previous 
time  step,  involved  into  the  non-linear  terms  (even 
after  linearization),  can  lead  to  instability. 

(iii)  Instability  due  to  the  conditions  at  the  interface. 
This  may  be  illustrated  by  the  simplest  example  of  Eq. 
(11a).  If  the  lengths  of  the  space  steps  inside  the 
separator  and  the  electrode  are  equal  to  each  other 
and  the  number  of  interface  points  is  n,  then  this 
condition  can  be  approximated  as  D*{f(Cn-Cn-1)  = 
Dsca(C„  +  1  —  C„).  From  that  expression  we  can  combine 
a  ‘production’  term  d2C/dz2  =  (Cn+1-2C„  +  Cn-1)/Ah2 
which,  roughly  speaking,  must  correspond  to  the  con¬ 
sumption  of  matter  in  the  interface  region.  However, 
it  does  not.  When  D*tt=D%n,  we  have  92C/9z2=0,  while 
the  production  must  vanish  under  completely  different 
conditions  —  the  absence  of  an  available  active  surface 
in  the  boundary  region  of  the  electrode.  Therefore, 
there  is  an  inaccuracy,  generally  inversely  proportional 
to  the  number  of  grid  points,  but  increasing  significantly 
in  typical  situations  when  a  discharge  mostly  takes  place 
in  the  electrode  in  this,  the  most  accessible,  part  of 
electrode.  Such  an  inaccuracy  of  10  or  even  20%  can 
cause  a  loss  of  stability. 

We  have  to  note  that  such  an  effect  is  not  directly 
correlated  with  the  accuracy  of  the  approximation  of 
the  boundary  conditions.  Moreover,  as  it  was  first  noted 
by  Arakawa  [20],  three-point  approximations  for  bound¬ 
ary  conditions,  as  a  rule,  have  much  worse  stability 
than  the  two-point  ones.  This  is  probably  why  the 
algorithm  for  solution  with  five-diagonal  matrix  de¬ 
scribed  in  Ref.  [17],  turned  out  not  to  work  for  our 
problem. 

In  conclusion,  ‘interface  conditions’  must  be  avoided 
whenever  it  is  possible. 

(iv)  Instability  due  to  boundary  conditions.  The  stability 
of  a  numerical  solution  also  depends  upon  the  ratio 
of  the  numbers  of  boundary  conditions  for  functions 
and  derivatives.  The  boundary  conditions  for  functions 
improve  stability,  while  boundary  conditions  for  deriv¬ 


atives  make  the  stability  worse.  It  is  connected  with 
some  fundamental  properties  of  the  original  Thomson 
algorithm  and  that  cannot  be  analysed  here  in  detail. 
Note,  however,  that  the  direction  in  which  the  differ¬ 
ential  equation  is  solved  is  also  important,  although 
from  a  mathematical  point  of  view  these  problems  are 
equivalent.  If  we  start  from  a  definite  value  of  the 
function  on  one  edge  and  meet  the  conditions  for 
derivative  on  the  other  end,  we  will  have  a  better  result 
than  in  the  opposite  case,  when  we  start  from  the 
derivative. 

In  our  problem  there  are  no  boundary  conditions 
for  porosity  and  all  boundary  conditions  for  concen¬ 
tration  and  potential  are  conditions  for  their  derivatives, 
which  leads  to  a  formal  instability  of  the  linearized 
homogeneous  system.  This  phenomenon  can  be  easily 
understood.  If  a  differential  equation  does  not  include 
the  function  itself  but  only  its  first  and  second  derivatives, 
a  function,  different  by  a  constant,  will  be  a  solution 
of  this  equation  and  boundary  conditions.  Moreover, 
two  boundary  conditions  for  derivatives  can  contradict 
each  other.  Therefore,  we  have  to  admit,  that  a  straight¬ 
forward  application  of  Newman’s  technique  [18]  to  the 
linearized  system  sometimes  can  involve  serious  math¬ 
ematical  and  numerical  difficulties.  Note  that,  strictly 
speaking,  our  speculation  is  valid  only  for  linearized 
elliptic  systems.  Of  course,  the  full  system  of  parabolic 
equations  (Eqs.  (1)  to  (8))  with  initial  conditions  (Eq. 
(9))  and  boundary  conditions  (Eqs.  (10)  to  (12))  are 
self-consistent  from  a  mathematical  point  of  view.  The 
values  of  function  C  and  X  cannot  be  set  arbitrarily, 
because  they  depend  on  the  initial  conditions  as  well. 

The  conservativity  of  a  numerical  scheme  is  a  measure 
of  how  long  it  is  able  to  describe  the  model  system 
with  necessary  accuracy.  In  some  sense  conservativity 
is  a  result  of  a  good  precision  and  stability.  At  the 
same  time  conservativity  is  a  wider  concept,  because 
it  requires  that  some  additional  parameters  keep  their 
values.  For  example,  on  one  hand,  in  our  case  of 
galvanostatic  discharge,  the  total  amounts  of  H2S04, 
Pb  and  PbOz  must  change  linearly  with  time,  because 
the  electrochemical  reactions  are  subject  to  Faraday’s 
law.  On  the  other  hand,  these  additional  conditions 
do  not  follow  directly  from  discretised  equations  and 
boundary  conditions. 

There  are  two  main  reasons  why  the  system  can  loose 
conservativity: 

(i)  Numerical.  All  the  above-mentioned  reasons  for 
non-stability  can  interfere  with  conservativity  as  well. 
We  also  have  to  note  that  there  is  no  general  algorithm 
for  the  construction  of  a  conservative  numerical  scheme 
for  non-linear  differential  equations,  as  those  appearing 
in  the  used  model.  Actually,  the  only  way  to  improve 
conservativity  is  to  perform  the  calculations  at  each 
time  level  with  sufficient  precision.  Our  experience 
shows  that  although  there  are  no  theoretical  restrictions 
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on  the  choice  of  a  grid  step  inside  each  region,  a 
uniform  space  step  for  all  regions  gives  a  much  more 
conservative  numerical  scheme.  Therefore,  in  the  com¬ 
puter  experiments  presented  in  this  paper,  we  always 
use  a  uniform  grid,  choosing  the  number  of  grid  points 
in  such  a  way  that  interface  points  placed  halfway 
between  two  grid  points  are  as  accurate  as  possible. 

(ii)  Model.  Some  model  approximations  also  can  lead 
to  non-conservative  solutions.  For  example,  it  is  gen¬ 
erally  accepted  that  the  transference  number  depends 
only  on  concentration  and  does  not  depend  on  current 
density.  Of  course,  this  dependence  is  rather  slow,  but 
integration  of  it  with  time,  at  a  high  gradient  of  current 
density,  can  give  us  a  significant  departure  for  con- 
servativity. 

Taking  into  account  all  the  above-mentioned  factors, 
we  have  constructed  our  own  numerical  algorithm.  This 
algorithm,  although  it  is  not  as  transparent  as  those 
of  Newman  [19]  and  White  [17],  gives  us  the  possibility 
of  a  stable,  accurate  and  fast  modelling  for  any  value 
of  the  input  parameters  up  to  a  high  degree  of  discharge. 
The  minimal  dimensionless  concentration  inside  the 
electrode  can  be  as  small  as  0.005. 
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